function llhess = HessPhat(theta,~,Y,X,Z,YY,V_R,V_US,V_UK,V_FRN,V_RUS,V_CHN,U0,intervEq,DMES,reb,mpow)

N=size(X,1);
KZ=size(Z,2)/5;
NS=size(V_R,2);
nBeta = length(reb{1});
nPhi1 = length(mpow{1});
nPhi2 = length(mpow{2});
nGamma1 = length(reb{2});
nGamma2 = length(reb{3});
nChi1 = length(mpow{3});
nChi2 = length(mpow{4});
nYY=size(YY,1);

Xrebel = X(:,reb{1});
Xmpow1 = X(:,mpow{1});
Xmpow2 = X(:,mpow{2});

ZUS = Z(:,1:KZ);
ZUK = Z(:,KZ+(1:KZ)); 
ZFRN = Z(:,2*KZ+(1:KZ));
ZRUS = Z(:,3*KZ+(1:KZ));
ZCHN = Z(:,4*KZ+(1:KZ));

ZrebelUS1 = ZUS(:,reb{2}); ZrebelUS2 = ZUS(:,reb{3});
ZrebelUK1 = ZUK(:,reb{2}); ZrebelUK2 = ZUK(:,reb{3});
ZrebelFRN1 = ZFRN(:,reb{2}); ZrebelFRN2 = ZFRN(:,reb{3});
ZrebelRUS1 = ZRUS(:,reb{2}); ZrebelRUS2 = ZRUS(:,reb{3});
ZrebelCHN1 = ZCHN(:,reb{2}); ZrebelCHN2 = ZCHN(:,reb{3});
ZmpowUS1 = ZUS(:,mpow{3}); ZmpowUS2 = ZUS(:,mpow{4});
ZmpowUK1 = ZUK(:,mpow{3}); ZmpowUK2 = ZUK(:,mpow{4});
ZmpowFRN1 = ZFRN(:,mpow{3}); ZmpowFRN2 = ZFRN(:,mpow{4});
ZmpowRUS1 = ZRUS(:,mpow{3}); ZmpowRUS2 = ZRUS(:,mpow{4});
ZmpowCHN1 = ZCHN(:,mpow{3}); ZmpowCHN2 = ZCHN(:,mpow{4});

beta=theta(1:nBeta,1);

phi_US1=theta(nBeta+[1,5+(1:nPhi1-1)],1);
phi_UK1=theta(nBeta+[2,5+(1:nPhi1-1)],1);
phi_FRN1=theta(nBeta+[3,5+(1:nPhi1-1)],1);
phi_RUS1=theta(nBeta+[4,5+(1:nPhi1-1)],1);
phi_CHN1=theta(nBeta+[5,5+(1:nPhi1-1)],1);
phi_US2=theta(nBeta+4+nPhi1+[1,5+(1:nPhi2-1)],1);
phi_UK2=theta(nBeta+4+nPhi1+[2,5+(1:nPhi2-1)],1);
phi_FRN2=theta(nBeta+4+nPhi1+[3,5+(1:nPhi2-1)],1);
phi_RUS2=theta(nBeta+4+nPhi1+[4,5+(1:nPhi2-1)],1);
phi_CHN2=theta(nBeta+4+nPhi1+[5,5+(1:nPhi2-1)],1);
    
gamma_US1=theta(nBeta+4+nPhi1+4+nPhi2+1,1);
gamma_UK1=theta(nBeta+4+nPhi1+4+nPhi2+2,1);
gamma_FRN1=theta(nBeta+4+nPhi1+4+nPhi2+3,1);
gamma_RUS1=theta(nBeta+4+nPhi1+4+nPhi2+4,1);
gamma_CHN1=theta(nBeta+4+nPhi1+4+nPhi2+5,1);
gamma01=theta(nBeta+4+nPhi1+4+nPhi2+5+(1:nGamma1),1);
gamma_US2=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+1,1);
gamma_UK2=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+2,1);
gamma_FRN2=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+3,1);
gamma_RUS2=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+4,1);
gamma_CHN2=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5,1);
gamma02=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+(1:nGamma2),1);

chi1=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+(1:nChi1),1);
chi2=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+(1:nChi2),1);

deltaS_US_UK=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+1,1);
deltaS_US_FRN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+2,1);
deltaS_US_RUS=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+3,1);
deltaS_US_CHN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+4,1);
deltaS_UK_FRN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+5,1);
deltaS_UK_RUS=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+6,1);
deltaS_UK_CHN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+7,1);
deltaS_FRN_RUS=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+8,1);
deltaS_FRN_CHN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+9,1);
deltaS_RUS_CHN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10,1);
deltaO_US_UK=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+1,1);
deltaO_US_FRN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+2,1);
deltaO_US_RUS=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+3,1);
deltaO_US_CHN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+4,1);
deltaO_UK_FRN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+5,1);
deltaO_UK_RUS=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+6,1);
deltaO_UK_CHN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+7,1);
deltaO_FRN_RUS=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+8,1);
deltaO_FRN_CHN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+9,1);
deltaO_RUS_CHN=theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+10+10,1);

lambda = theta(nBeta+4+nPhi1+4+nPhi2+5+nGamma1+5+nGamma2+nChi1+nChi2+21:end,1);

llhess=zeros(length(theta));

for n=1:N
    
    u_R = YY(:,1) .* (Xrebel(n,:)*beta + ...
        (YY(:,2)==1)*(gamma_US1+ZrebelUS1(n,:)*gamma01) + ...
        (YY(:,2)==2)*(gamma_US2+ZrebelUS2(n,:)*gamma02) + ...
        (YY(:,3)==1)*(gamma_UK1+ZrebelUK1(n,:)*gamma01) + ...
        (YY(:,3)==2)*(gamma_UK2+ZrebelUK2(n,:)*gamma02) + ...
        (YY(:,4)==1)*(gamma_FRN1+ZrebelFRN1(n,:)*gamma01) + ...
        (YY(:,4)==2)*(gamma_FRN2+ZrebelFRN2(n,:)*gamma02) + ...
        (YY(:,5)==1)*(gamma_RUS1+ZrebelRUS1(n,:)*gamma01) + ...
        (YY(:,5)==2)*(gamma_RUS2+ZrebelRUS2(n,:)*gamma02) + ...
        (YY(:,6)==1)*(gamma_CHN1+ZrebelCHN1(n,:)*gamma01) + ...
        (YY(:,6)==2)*(gamma_CHN2+ZrebelCHN2(n,:)*gamma02));
    
    u_US = (YY(:,2)==1) .* (Xmpow1(n,:)*phi_US1 + ZmpowUS1(n,:)*chi1 + ...
        (YY(:,3)==1)*deltaS_US_UK + (YY(:,3)==2)*deltaO_US_UK + ...
        (YY(:,4)==1)*deltaS_US_FRN + (YY(:,4)==2)*deltaO_US_FRN + ...
        (YY(:,5)==1)*deltaS_US_RUS + (YY(:,5)==2)*deltaO_US_RUS + ...
        (YY(:,6)==1)*deltaS_US_CHN + (YY(:,6)==2)*deltaO_US_CHN) + ...
        (YY(:,2)==2) .* (Xmpow2(n,:)*phi_US2 + ZmpowUS2(n,:)*chi2 + ...
        (YY(:,3)==1)*deltaO_US_UK + (YY(:,3)==2)*deltaS_US_UK + ...
        (YY(:,4)==1)*deltaO_US_FRN + (YY(:,4)==2)*deltaS_US_FRN + ...
        (YY(:,5)==1)*deltaO_US_RUS + (YY(:,5)==2)*deltaS_US_RUS + ...
        (YY(:,6)==1)*deltaO_US_CHN + (YY(:,6)==2)*deltaS_US_CHN);
    
    u_UK = (YY(:,3)==1) .* (Xmpow1(n,:)*phi_UK1 + ZmpowUK1(n,:)*chi1 + ...
        (YY(:,2)==1)*deltaS_US_UK + (YY(:,2)==2)*deltaO_US_UK + ...
        (YY(:,4)==1)*deltaS_UK_FRN + (YY(:,4)==2)*deltaO_UK_FRN + ...
        (YY(:,5)==1)*deltaS_UK_RUS + (YY(:,5)==2)*deltaO_UK_RUS + ...
        (YY(:,6)==1)*deltaS_UK_CHN + (YY(:,6)==2)*deltaO_UK_CHN) + ...
        (YY(:,3)==2) .* (Xmpow2(n,:)*phi_UK2 + ZmpowUK2(n,:)*chi2 + ...
        (YY(:,2)==1)*deltaO_US_UK + (YY(:,2)==2)*deltaS_US_UK + ...
        (YY(:,4)==1)*deltaO_UK_FRN + (YY(:,4)==2)*deltaS_UK_FRN + ...
        (YY(:,5)==1)*deltaO_UK_RUS + (YY(:,5)==2)*deltaS_UK_RUS + ...
        (YY(:,6)==1)*deltaO_UK_CHN + (YY(:,6)==2)*deltaS_UK_CHN);

    u_FRN = (YY(:,4)==1) .* (Xmpow1(n,:)*phi_FRN1 + ZmpowFRN1(n,:)*chi1 + ...
        (YY(:,2)==1)*deltaS_US_FRN + (YY(:,2)==2)*deltaO_US_FRN + ...
        (YY(:,3)==1)*deltaS_UK_FRN + (YY(:,3)==2)*deltaO_UK_FRN + ...
        (YY(:,5)==1)*deltaS_FRN_RUS + (YY(:,5)==2)*deltaO_FRN_RUS + ...
        (YY(:,6)==1)*deltaS_FRN_CHN + (YY(:,6)==2)*deltaO_FRN_CHN) + ...
        (YY(:,4)==2) .* (Xmpow2(n,:)*phi_FRN2 + ZmpowFRN2(n,:)*chi2 + ...
        (YY(:,2)==1)*deltaO_US_FRN + (YY(:,2)==2)*deltaS_US_FRN + ...
        (YY(:,3)==1)*deltaO_UK_FRN + (YY(:,3)==2)*deltaS_UK_FRN + ...
        (YY(:,5)==1)*deltaO_FRN_RUS + (YY(:,5)==2)*deltaS_FRN_RUS + ...
        (YY(:,6)==1)*deltaO_FRN_CHN + (YY(:,6)==2)*deltaS_FRN_CHN);
    
    u_RUS = (YY(:,5)==1) .* (Xmpow1(n,:)*phi_RUS1 + ZmpowRUS1(n,:)*chi1 + ...
        (YY(:,2)==1)*deltaS_US_RUS + (YY(:,2)==2)*deltaO_US_RUS + ...
        (YY(:,3)==1)*deltaS_UK_RUS + (YY(:,3)==2)*deltaO_UK_RUS + ...
        (YY(:,4)==1)*deltaS_FRN_RUS + (YY(:,4)==2)*deltaO_FRN_RUS + ...
        (YY(:,6)==1)*deltaS_RUS_CHN + (YY(:,6)==2)*deltaO_RUS_CHN) + ...
        (YY(:,5)==2) .* (Xmpow2(n,:)*phi_RUS2 + ZmpowRUS2(n,:)*chi2 + ...
        (YY(:,2)==1)*deltaO_US_RUS + (YY(:,2)==2)*deltaS_US_RUS + ...
        (YY(:,3)==1)*deltaO_UK_RUS + (YY(:,3)==2)*deltaS_UK_RUS + ...
        (YY(:,4)==1)*deltaO_FRN_RUS + (YY(:,4)==2)*deltaS_FRN_RUS + ...
        (YY(:,6)==1)*deltaO_RUS_CHN + (YY(:,6)==2)*deltaS_RUS_CHN);

    u_CHN = (YY(:,6)==1) .* (Xmpow1(n,:)*phi_CHN1 + ZmpowCHN1(n,:)*chi1 + ...
        (YY(:,2)==1)*deltaS_US_CHN + (YY(:,2)==2)*deltaO_US_CHN + ...
        (YY(:,3)==1)*deltaS_UK_CHN + (YY(:,3)==2)*deltaO_UK_CHN + ...
        (YY(:,4)==1)*deltaS_FRN_CHN + (YY(:,4)==2)*deltaO_FRN_CHN + ...
        (YY(:,5)==1)*deltaS_RUS_CHN + (YY(:,5)==2)*deltaO_RUS_CHN) + ...
        (YY(:,6)==2) .* (Xmpow2(n,:)*phi_CHN2 + ZmpowCHN2(n,:)*chi2 + ...
        (YY(:,2)==1)*deltaO_US_CHN + (YY(:,2)==2)*deltaS_US_CHN + ...
        (YY(:,3)==1)*deltaO_UK_CHN + (YY(:,3)==2)*deltaS_UK_CHN + ...
        (YY(:,4)==1)*deltaO_FRN_CHN + (YY(:,4)==2)*deltaS_FRN_CHN + ...
        (YY(:,5)==1)*deltaO_RUS_CHN + (YY(:,5)==2)*deltaS_RUS_CHN);
       
    V = [V_R(:,:,n);V_US(:,:,n);V_UK(:,:,n);V_FRN(:,:,n);V_RUS(:,:,n);V_CHN(:,:,n)]';
    U = [u_R;u_US;u_UK;u_FRN;u_RUS;u_CHN]';
    U0n = reshape(U0(:,:,n),6*nYY,1)';
    weights = exp(-sum((V-U).^2,2)/2 + sum((V-U0n).^2,2)/2) / NS;    
    
    % Gradient w.r.t. payoffs
    Du = zeros(6*nYY,length(theta)-length(lambda));
    % rebel payoffs  
    Du(1:nYY,[1:nBeta,nBeta+4+nPhi1+4+nPhi2+(1:2*5+nGamma1+nGamma2)]) = ...
    	[YY(:,1)*Xrebel(n,:), YY(:,2:end)==1, (YY(:,2)==1)*ZrebelUS1(n,:) + ...
        (YY(:,3)==1)*ZrebelUK1(n,:) + (YY(:,4)==1)*ZrebelFRN1(n,:) + ...
        (YY(:,5)==1)*ZrebelRUS1(n,:) + (YY(:,6)==1)*ZrebelCHN1(n,:), ...
        YY(:,2:end)==2, (YY(:,2)==2)*ZrebelUS2(n,:) + ...
        (YY(:,3)==2)*ZrebelUK2(n,:) + (YY(:,4)==2)*ZrebelFRN2(n,:) + ...
        (YY(:,5)==2)*ZrebelRUS2(n,:) + (YY(:,6)==2)*ZrebelCHN2(n,:)];
    % US payoffs
    Du(nYY+(1:nYY),[nBeta+[1,5+(1:nPhi1-1)],nBeta+4+nPhi1+[1,5+(1:nPhi2-1)]...
    	nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+(1:nChi1+nChi2),...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+(1:4),...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+10+(1:4)]) = ...
        (YY(:,2)==1) .* [YY(:,1)*[Xmpow1(n,:), zeros(1,nPhi2), ZmpowUS1(n,:), zeros(1,nChi2)], ...
        YY(:,3:6)==1, YY(:,3:6)==2] + ...
        (YY(:,2)==2) .* [YY(:,1)*[zeros(1,nPhi1), Xmpow2(n,:), zeros(1,nChi1), ZmpowUS2(n,:)], ...
        YY(:,3:6)==2, YY(:,3:6)==1];
	% UK
    Du(2*nYY+(1:nYY),[nBeta+[2,5+(1:nPhi1-1)],nBeta+4+nPhi1+[2,5+(1:nPhi2-1)]...
    	nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+(1:nChi1+nChi2),...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+[1,5:7],...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+10+[1,5:7]]) = ...
        (YY(:,3)==1) .* [YY(:,1)*[Xmpow1(n,:), zeros(1,nPhi2), ZmpowUK1(n,:), zeros(1,nChi2)], ...
        YY(:,[2,4:6])==1, YY(:,[2,4:6])==2] + ...
        (YY(:,3)==2) .* [YY(:,1)*[zeros(1,nPhi1), Xmpow2(n,:), zeros(1,nChi1), ZmpowUK2(n,:)], ...
        YY(:,[2,4:6])==2, YY(:,[2,4:6])==1];
	% FRN 
    Du(3*nYY+(1:nYY),[nBeta+[3,5+(1:nPhi1-1)],nBeta+4+nPhi1+[3,5+(1:nPhi2-1)]...
    	nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+(1:nChi1+nChi2),...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+[2,5,8:9],...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+10+[2,5,8:9]]) = ...
        (YY(:,4)==1) .* [YY(:,1)*[Xmpow1(n,:), zeros(1,nPhi2), ZmpowFRN1(n,:), zeros(1,nChi2)], ...
        YY(:,[2:3,5:6])==1, YY(:,[2:3,5:6])==2] + ...
        (YY(:,4)==2) .* [YY(:,1)*[zeros(1,nPhi1), Xmpow2(n,:), zeros(1,nChi1), ZmpowFRN2(n,:)], ...
        YY(:,[2:3,5:6])==2, YY(:,[2:3,5:6])==1];
	% RUS
    Du(4*nYY+(1:nYY),[nBeta+[4,5+(1:nPhi1-1)],nBeta+4+nPhi1+[4,5+(1:nPhi2-1)]...
    	nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+(1:nChi1+nChi2),...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+[3,6,8,10],...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+10+[3,6,8,10]]) = ...
        (YY(:,5)==1) .* [YY(:,1)*[Xmpow1(n,:), zeros(1,nPhi2), ZmpowRUS1(n,:), zeros(1,nChi2)], ...
        YY(:,[2:4,6])==1, YY(:,[2:4,6])==2] + ...
        (YY(:,5)==2) .* [YY(:,1)*[zeros(1,nPhi1), Xmpow2(n,:), zeros(1,nChi1), ZmpowRUS2(n,:)], ...
        YY(:,[2:4,6])==2, YY(:,[2:4,6])==1];
	% CHN
    Du(5*nYY+(1:nYY),[nBeta+[5,5+(1:nPhi1-1)],nBeta+4+nPhi1+[5,5+(1:nPhi2-1)]...
    	nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+(1:nChi1+nChi2),...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+[4,7,9:10],...
        nBeta+4+nPhi1+4+nPhi2+2*5+nGamma1+nGamma2+nChi1+nChi2+10+[4,7,9:10]]) = ...
        (YY(:,6)==1) .* [YY(:,1)*[Xmpow1(n,:), zeros(1,nPhi2), ZmpowCHN1(n,:), zeros(1,nChi2)], ...
        YY(:,2:5)==1, YY(:,2:5)==2] + ...
        (YY(:,6)==2) .* [YY(:,1)*[zeros(1,nPhi1), Xmpow2(n,:), zeros(1,nChi1), ZmpowCHN2(n,:)], ...
        YY(:,2:5)==2, YY(:,2:5)==1];

    Pra = zeros(NS, 1);
    Grad = zeros(NS, length(lambda));
    H22=zeros(length(lambda));
	for b=1:NS
    	[Pra(b,1), Grad(b,:), h22] = profprob(Y(n,1), DMES{b,n}, intervEq{b,n}, V_R(1,b,n), lambda);
        H22 = H22 + weights(b,1) * h22;
	end
    
    ln = weights' * Pra;
    lg1 = weights' * (Pra .* (V-U) * Du) / ln;
    lg2 = weights' * Grad / ln;
    
    H11 = - Du' * Du + (weights .* (V-U) * Du)' * (Pra .* (V-U) * Du) / ln - lg1' * lg1;
    H12 = (weights .* (V-U) * Du)' * Grad / ln - lg1' * lg2;
    H22 = H22 / ln - lg2' * lg2;
    
    llhess = llhess - [H11, H12; H12', H22];
        
end         
    


